//
// Function mv1Eval() to compute the MV1 b-tagging weight of a jet.
// Signature: double mv1Eval(double w_IP3D, double w_SV1, double w_JetFitterCombNN, 
//                           double jet_pt, double jet_eta)
//
// Inputs are: IP3D tagging weight, SV1 tagging weight, JetFitterCombNN tagging weight, 
//             jet pT (in MeV) and jet eta.
//
// To use it in a ROOT macro, compile it: gROOT->ProcessLine(".L MV1.C+");
//                            or include it: #include "MV1.C" 
//
// 2011/11/19 L.Vacavant
//
//----------------------------------------------------------------------------------------------

// Class: ReadMLP_ANN
// Automatically generated by MethodBase::MakeClass
//

/* configuration options =====================================================

#GEN -*-*-*-*-*-*-*-*-*-*-*- general info -*-*-*-*-*-*-*-*-*-*-*-

Method         : MLP::MLP_ANN
TMVA Release   : 4.1.2         [262402]
ROOT Release   : 5.30/03       [335363]
Creator        : vacavant
Date           : Wed Nov 16 19:23:52 2011
Host           : Linux lxbuild148.cern.ch 2.6.18-194.26.1.el5 #1 SMP Wed Nov 10 09:45:46 CET 2010 x86_64 x86_64 x86_64 GNU/Linux
Dir            : /atlas/vacavant/work/btagPerfR17/micront_analyse/ANNweighted/wptcat_nosmt_jfc_mlp32
Training events: 5165297
Analysis type  : [Classification]


#OPT -*-*-*-*-*-*-*-*-*-*-*-*- options -*-*-*-*-*-*-*-*-*-*-*-*-

# Set by User:
NCycles: "3000" [Number of training cycles]
HiddenLayers: "3,2" [Specification of hidden layer architecture]
V: "True" [Verbose output (short form of "VerbosityLevel" below - overrides the latter one)]
H: "True" [Print method-specific help message]
TestRate: "5" [Test for overtraining performed at each #th epochs]
# Default:
NeuronType: "sigmoid" [Neuron activation function type]
RandomSeed: "1" [Random seed for initial synapse weights (0 means unique seed for each run; default value '1')]
EstimatorType: "MSE" [MSE (Mean Square Estimator) for Gaussian Likelihood or CE(Cross-Entropy) for Bernoulli Likelihood]
NeuronInputType: "sum" [Neuron input function type]
VerbosityLevel: "Verbose" [Verbosity level]
VarTransform: "None" [List of variable transformations performed before training, e.g., "D_Background,P_Signal,G,N_AllClasses" for: "Decorrelation, PCA-transformation, Gaussianisation, Normalisation, each for the given class of events ('AllClasses' denotes all events of all classes, if no class indication is given, 'All' is assumed)"]
CreateMVAPdfs: "False" [Create PDFs for classifier outputs (signal and background)]
IgnoreNegWeightsInTraining: "False" [Events with negative weights are ignored in the training (but are included for testing and performance evaluation)]
TrainingMethod: "BP" [Train with Back-Propagation (BP), BFGS Algorithm (BFGS), or Genetic Algorithm (GA - slower and worse)]
LearningRate: "2.000000e-02" [ANN learning rate parameter]
DecayRate: "1.000000e-02" [Decay rate for learning parameter]
EpochMonitoring: "False" [Provide epoch-wise monitoring plots according to TestRate (caution: causes big ROOT output file!)]
Sampling: "1.000000e+00" [Only 'Sampling' (randomly selected) events are trained each epoch]
SamplingEpoch: "1.000000e+00" [Sampling is used for the first 'SamplingEpoch' epochs, afterwards, all events are taken for training]
SamplingImportance: "1.000000e+00" [ The sampling weights of events in epochs which successful (worse estimator than before) are multiplied with SamplingImportance, else they are divided.]
SamplingTraining: "True" [The training sample is sampled]
SamplingTesting: "False" [The testing sample is sampled]
ResetStep: "50" [How often BFGS should reset history]
Tau: "3.000000e+00" [LineSearch "size step"]
BPMode: "sequential" [Back-propagation learning mode: sequential or batch]
BatchSize: "-1" [Batch size: number of events/batch, only set if in Batch Mode, -1 for BatchSize=number_of_events]
ConvergenceImprove: "1.000000e-30" [Minimum improvement which counts as improvement (<0 means automatic convergence check is turned off)]
ConvergenceTests: "-1" [Number of steps (without improvement) required for convergence (<0 means automatic convergence check is turned off)]
UseRegulator: "False" [Use regulator to avoid over-training]
UpdateLimit: "10000" [Maximum times of regulator update]
CalculateErrors: "False" [Calculates inverse Hessian matrix at the end of the training to be able to calculate the uncertainties of an MVA value]
WeightRange: "1.000000e+00" [Take the events for the estimator calculations from small deviations from the desired value to large deviations only over the weight range]
##


#VAR -*-*-*-*-*-*-*-*-*-*-*-* variables *-*-*-*-*-*-*-*-*-*-*-*-

NVar 4
ip3                           ip3                           ip3                           ip3                                                             'F'    [-30,100]
sv1                           sv1                           sv1                           sv1                                                             'F'    [-30,11.1977300644]
jfc                           jfc                           jfc                           jfc                                                             'F'    [-8.91371440887,10.1707715988]
categ(pt,eta)                 categ_pt,eta_                 categ(pt,eta)                 categ(pt,eta)                                                   'I'    [13,58]
NSpec 8
smtfilt(smt)                  smtfilt_smt_                  smtfilt(smt)                  smtfilt(smt)                                                    'F'    [-5,12.7354049683]
evt                           evt                           evt                           I                                                               'F'    [1,1400000]
pt                            pt                            pt                            pt                                                              'F'    [20000.0078125,1311420.875]
eta                           eta                           eta                           eta                                                             'F'    [-2.49999952316,2.49999547005]
phi                           phi                           phi                           phi                                                             'F'    [-3.14159250259,3.1415913105]
cmb                           cmb                           cmb                           cmb                                                             'F'    [-21.9221591949,48.7108039856]
jft                           jft                           jft                           jft                                                             'F'    [-8.81255912781,5.88171577454]
label                         label                         label                         I                                                               'F'    [0,5]


============================================================================ */

#include <vector>
#include <cmath>
#include <string>
#include <iostream>

#ifndef IClassifierReader__def
#define IClassifierReader__def

class IClassifierReader {

 public:

   // constructor
   IClassifierReader() : fStatusIsClean( true ) {}
   virtual ~IClassifierReader() {}

   // return classifier response
   virtual double GetMvaValue( const std::vector<double>& inputValues ) const = 0;

   // returns classifier status
   bool IsStatusClean() const { return fStatusIsClean; }

 protected:

   bool fStatusIsClean;
};

#endif

class ReadMLP_ANN : public IClassifierReader {

 public:

   // constructor
   ReadMLP_ANN( std::vector<std::string>& theInputVars ) 
      : IClassifierReader(),
        fClassName( "ReadMLP_ANN" ),
        fNvars( 4 ),
        fIsNormalised( false )
   {      
      // the training input variables
      const char* inputVars[] = { "ip3", "sv1", "jfc", "categ(pt,eta)" };

      // sanity checks
      if (theInputVars.size() <= 0) {
         std::cout << "Problem in class \"" << fClassName << "\": empty input vector" << std::endl;
         fStatusIsClean = false;
      }

      if (theInputVars.size() != fNvars) {
         std::cout << "Problem in class \"" << fClassName << "\": mismatch in number of input values: "
                   << theInputVars.size() << " != " << fNvars << std::endl;
         fStatusIsClean = false;
      }

      // validate input variables
      for (size_t ivar = 0; ivar < theInputVars.size(); ivar++) {
         if (theInputVars[ivar] != inputVars[ivar]) {
            std::cout << "Problem in class \"" << fClassName << "\": mismatch in input variable names" << std::endl
                      << " for variable [" << ivar << "]: " << theInputVars[ivar].c_str() << " != " << inputVars[ivar] << std::endl;
            fStatusIsClean = false;
         }
      }

      // initialize min and max vectors (for normalisation)
      fVmin[0] = -30;
      fVmax[0] = 100;
      fVmin[1] = -30;
      fVmax[1] = 11.1977300643921;
      fVmin[2] = -8.91371440887451;
      fVmax[2] = 10.1707715988159;
      fVmin[3] = 13;
      fVmax[3] = 58;

      // initialize input variable types
      fType[0] = 'F';
      fType[1] = 'F';
      fType[2] = 'F';
      fType[3] = 'I';

      // initialize constants
      Initialize();

   }

   // destructor
   virtual ~ReadMLP_ANN() {
      Clear(); // method-specific
   }

   // the classifier response
   // "inputValues" is a vector of input values in the same order as the 
   // variables given to the constructor
   double GetMvaValue( const std::vector<double>& inputValues ) const;

 private:

   // method-specific destructor
   void Clear();

   // common member variables
   const char* fClassName;

   const size_t fNvars;
   size_t GetNvar()           const { return fNvars; }
   char   GetType( int ivar ) const { return fType[ivar]; }

   // normalisation of input variables
   const bool fIsNormalised;
   bool IsNormalised() const { return fIsNormalised; }
   double fVmin[4];
   double fVmax[4];
   double NormVariable( double x, double xmin, double xmax ) const {
      // normalise to output range: [-1, 1]
      return 2*(x - xmin)/(xmax - xmin) - 1.0;
   }

   // type of input variable: 'F' or 'I'
   char   fType[4];

   // initialize internal variables
   void Initialize();
   double GetMvaValue__( const std::vector<double>& inputValues ) const;

   // private members (method specific)

   double ActivationFnc(double x) const;
   double OutputActivationFnc(double x) const;

   int fLayers;
   int fLayerSize[4];
   double fWeightMatrix0to1[4][5];   // weight matrix from layer 0 to 1
   double fWeightMatrix1to2[3][4];   // weight matrix from layer 1 to 2
   double fWeightMatrix2to3[1][3];   // weight matrix from layer 2 to 3

   double * fWeights[4];
};

inline void ReadMLP_ANN::Initialize()
{
   // build network structure
   fLayers = 4;
   fLayerSize[0] = 5; fWeights[0] = new double[5]; 
   fLayerSize[1] = 4; fWeights[1] = new double[4]; 
   fLayerSize[2] = 3; fWeights[2] = new double[3]; 
   fLayerSize[3] = 1; fWeights[3] = new double[1]; 
   // weight matrix from layer 0 to 1
   fWeightMatrix0to1[0][0] = -0.0645996654215314;
   fWeightMatrix0to1[1][0] = -0.00907307647258392;
   fWeightMatrix0to1[2][0] = 0.777218896879578;
   fWeightMatrix0to1[0][1] = 3.349292717926;
   fWeightMatrix0to1[1][1] = 0.381815537920455;
   fWeightMatrix0to1[2][1] = -1.43608372700896;
   fWeightMatrix0to1[0][2] = 4.20275878680485;
   fWeightMatrix0to1[1][2] = 1.03060385959224;
   fWeightMatrix0to1[2][2] = -1.27341173103263;
   fWeightMatrix0to1[0][3] = -0.057349339802193;
   fWeightMatrix0to1[1][3] = -0.0256639296849042;
   fWeightMatrix0to1[2][3] = -1.214375825434;
   fWeightMatrix0to1[0][4] = 9.14595969016655;
   fWeightMatrix0to1[1][4] = 0.715588833083844;
   fWeightMatrix0to1[2][4] = -0.677603487608217;
   // weight matrix from layer 1 to 2
   fWeightMatrix1to2[0][0] = 11.3008892808662;
   fWeightMatrix1to2[1][0] = -18.8208466881222;
   fWeightMatrix1to2[0][1] = 2.3515898891914;
   fWeightMatrix1to2[1][1] = -2.07980454105222;
   fWeightMatrix1to2[0][2] = 0.0667825242021229;
   fWeightMatrix1to2[1][2] = -0.579628050791058;
   fWeightMatrix1to2[0][3] = -14.4738288723496;
   fWeightMatrix1to2[1][3] = -3.54581799346469;
   // weight matrix from layer 2 to 3
   fWeightMatrix2to3[0][0] = 2.73279555773931;
   fWeightMatrix2to3[0][1] = -4.97404654862053;
   fWeightMatrix2to3[0][2] = 0.162561097335823;
}

inline double ReadMLP_ANN::GetMvaValue__( const std::vector<double>& inputValues ) const
{
   if (inputValues.size() != (unsigned int)fLayerSize[0]-1) {
      std::cout << "Input vector needs to be of size " << fLayerSize[0]-1 << std::endl;
      return 0;
   }

   for (int l=0; l<fLayers; l++)
      for (int i=0; i<fLayerSize[l]; i++) fWeights[l][i]=0;

   for (int l=0; l<fLayers-1; l++)
      fWeights[l][fLayerSize[l]-1]=1;

   for (int i=0; i<fLayerSize[0]-1; i++)
      fWeights[0][i]=inputValues[i];

   // layer 0 to 1
   for (int o=0; o<fLayerSize[1]-1; o++) {
      for (int i=0; i<fLayerSize[0]; i++) {
         double inputVal = fWeightMatrix0to1[o][i] * fWeights[0][i];
         fWeights[1][o] += inputVal;
      }
      fWeights[1][o] = ActivationFnc(fWeights[1][o]);
   }
   // layer 1 to 2
   for (int o=0; o<fLayerSize[2]-1; o++) {
      for (int i=0; i<fLayerSize[1]; i++) {
         double inputVal = fWeightMatrix1to2[o][i] * fWeights[1][i];
         fWeights[2][o] += inputVal;
      }
      fWeights[2][o] = ActivationFnc(fWeights[2][o]);
   }
   // layer 2 to 3
   for (int o=0; o<fLayerSize[3]; o++) {
      for (int i=0; i<fLayerSize[2]; i++) {
         double inputVal = fWeightMatrix2to3[o][i] * fWeights[2][i];
         fWeights[3][o] += inputVal;
      }
      fWeights[3][o] = OutputActivationFnc(fWeights[3][o]);
   }

   return fWeights[3][0];
}

double ReadMLP_ANN::ActivationFnc(double x) const {
   // sigmoid
   return 1.0/(1.0+exp(-x));
}
double ReadMLP_ANN::OutputActivationFnc(double x) const {
   // identity
   return x;
}
   
// Clean up
inline void ReadMLP_ANN::Clear() {
}

inline double ReadMLP_ANN::GetMvaValue( const std::vector<double>& inputValues ) const {
      // classifier response value
      double retval = 0;

      // classifier response, sanity check first
      if (!IsStatusClean()) {
         std::cout << "Problem in class \"" << fClassName << "\": cannot return classifier response"
                   << " because status is dirty" << std::endl;
         retval = 0;
      }
      else {
         if (IsNormalised()) {
            // normalise variables
            std::vector<double> iV;
            int ivar = 0;
            for (std::vector<double>::const_iterator varIt = inputValues.begin();
                 varIt != inputValues.end(); varIt++, ivar++) {
               iV.push_back(NormVariable( *varIt, fVmin[ivar], fVmax[ivar] ));
            }
            retval = GetMvaValue__( iV );
         }
         else {
            retval = GetMvaValue__( inputValues );
         }
      }

      return retval;
}

//----------------------------------------------------------------------------------------------
// User part:
#include <algorithm>

IClassifierReader* _mlpResponse;

// helper functions to define the jet category: 
static const int njptbin = 11;
static const double _jptbins[njptbin] = { 15., 30., 45., 60., 100., 140., 180., 220., 260., 310., 500. };
std::vector<double> _vptbins(_jptbins,_jptbins+njptbin);
static const int netabin = 5;
static const double _etabins[netabin] = { 0., 0.6, 1.2, 1.8, 2.5 };
std::vector<double> _vetabins(_etabins,_etabins+netabin);
int findptbin(double pt) {
  std::vector<double>::iterator pos = std::lower_bound(_vptbins.begin(), 
                                                       _vptbins.end(), pt);
  return (pos-_vptbins.begin());
}
int findetabin(double eta) {
  std::vector<double>::iterator pos = std::lower_bound(_vetabins.begin(), 
                                                       _vetabins.end(), eta);
  return (pos-_vetabins.begin());
}
int categ(double pt, double eta) {
  int cat = 0;
  double npt = pt/1000.; // input in MeV, internally in GeV
  double neta = fabs(eta)+0.00000001; // to get same results for boundaries as Root FindBin()
  int binx = findptbin(npt);
  int biny = findetabin(neta);
  cat = binx + (njptbin+1)*biny;
  return cat;
}

// initialization of the network response function: do not call it yourself
void mv1Init() {
    std::vector<std::string> inputVars;
    inputVars.push_back( "ip3" );
    inputVars.push_back( "sv1" );
    inputVars.push_back( "jfc" );
    inputVars.push_back( "categ(pt,eta)" );
    _mlpResponse = new ReadMLP_ANN( inputVars );
}

// the only function to be called by the user:
double mv1Eval(double w_IP3D, double w_SV1, double w_JetFitterCombNN, double jet_pt, double jet_eta) {
  if(_mlpResponse==0) mv1Init();
  double mv1 = 0.;
  if(_mlpResponse) {
    std::vector<double> inputVec(4);
    inputVec[0] = w_IP3D;
    inputVec[1] = w_SV1;
    inputVec[2] = w_JetFitterCombNN;
    inputVec[3] = (double)(categ(jet_pt,jet_eta));
    mv1 = _mlpResponse->GetMvaValue( inputVec );
  }
  return mv1;
}

